# [dataset]
data_files = ["./datasets/BikeNYC/NYC14_M16x8_T60_NewEnd.h5"]
holiday_file = ""
meteorol_file = ""
T = 24  # time step of one hour
len_closeness = 3
len_period = 1
len_trend = 1
period_interval = 1
trend_interval = 7
len_test = 672
use_meta = True

# [model]
nb_flow = 2
map_height = 8
map_width = 16
nb_residual_unit = 4

# [learning]
epochs = 50
batch_size = 32
learning_rate = 0.0002

import tables
import numpy as np

def remove_incomplete_days(dat: tables.file.File, T: int) -> tuple[np.ndarray, np.ndarray]:
    # 20140425 has 24 timestamps, which does not appear in `incomplete_days` in the original implementation.
    # So I reimplemented it in a different way.
    data = dat.root.data.read()
    timestamps = dat.root.date.read().astype(str)

    dates, values = np.vstack(
        np.frompyfunc(lambda x: (x[:8], x[8:]), 1, 2)(timestamps)
    )
    # label encoding
    uniq_dates, labels = np.unique(dates, return_inverse=True)
    # groupby("labels")["values"].sum() != sum(range(1, 49))
    incomplete_days = uniq_dates[
        np.where(np.bincount(labels, values.astype(int)) != sum(range(1, (T + 1))))[
            0
        ]
    ]
    del_idx = np.where(np.isin(dates, incomplete_days))[0]
    new_data = np.delete(data, del_idx, axis=0)
    new_timestamps = np.delete(timestamps, del_idx)
    return new_data, new_timestamps

dat=tables.open_file(data_files[0],mode="r")
data=dat.root.data.read()
data

array([[[[ 0.,  0.,  1., ...,  0.,  0.,  0.],
         [ 0.,  1.,  1., ...,  1.,  0.,  0.],
         [ 0.,  0.,  2., ...,  2.,  2.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  3.,  1.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  1.]],

        [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  0.,  1.,  0.],
         [ 0.,  0.,  0., ...,  2.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.]]],


       [[[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  0.,  1.,  2.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  1.]],

        [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  1., ...,  0.,  0.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  1.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  1.]]],


       [[[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  0.,  2.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.]],

        [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         [ 0.,  0.,  0., ...,  2.,  1.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  0.]]],


       ...,


       [[[ 0.,  0.,  5., ...,  0.,  0.,  0.],
         [ 0.,  2.,  6., ...,  6.,  0.,  0.],
         [ 0.,  5., 22., ..., 14.,  8.,  0.],
         ...,
         [ 0.,  0.,  4., ...,  9.,  5.,  2.],
         [ 0.,  0.,  3., ...,  7.,  3.,  1.],
         [ 0.,  0.,  0., ...,  4.,  6.,  8.]],

        [[ 0.,  0.,  2., ...,  0.,  0.,  0.],
         [ 0.,  2.,  2., ...,  9.,  0.,  0.],
         [ 0.,  0., 10., ..., 13.,  7.,  0.],
         ...,
         [ 0.,  0.,  3., ...,  7.,  4.,  1.],
         [ 0.,  0.,  2., ...,  3.,  2.,  2.],
         [ 0.,  0.,  0., ...,  2.,  1.,  1.]]],


       [[[ 0.,  0.,  5., ...,  0.,  0.,  0.],
         [ 0.,  1.,  4., ...,  5.,  0.,  0.],
         [ 0.,  1.,  7., ...,  3.,  7.,  0.],
         ...,
         [ 0.,  0.,  3., ...,  8.,  1.,  0.],
         [ 0.,  0.,  2., ...,  4.,  0.,  1.],
         [ 0.,  0.,  0., ...,  2.,  4.,  3.]],

        [[ 0.,  0.,  2., ...,  0.,  0.,  0.],
         [ 0.,  1.,  5., ...,  7.,  0.,  0.],
         [ 0.,  0.,  8., ...,  3.,  5.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  2.,  1.,  1.],
         [ 0.,  0.,  2., ...,  0.,  1.,  0.],
         [ 0.,  0.,  0., ...,  0.,  0.,  1.]]],


       [[[ 0.,  0.,  1., ...,  0.,  0.,  0.],
         [ 0.,  0.,  3., ...,  4.,  0.,  0.],
         [ 0.,  2.,  2., ...,  4.,  6.,  0.],
         ...,
         [ 0.,  0.,  0., ..., 12.,  0.,  2.],
         [ 0.,  0.,  1., ...,  3.,  1.,  0.],
         [ 0.,  0.,  0., ...,  2.,  2.,  2.]],

        [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
         [ 0.,  0.,  2., ...,  3.,  0.,  0.],
         [ 0.,  0.,  5., ...,  0.,  2.,  0.],
         ...,
         [ 0.,  0.,  0., ...,  1.,  0.,  0.],
         [ 0.,  0.,  1., ...,  1.,  0.,  0.],
         [ 0.,  0.,  0., ...,  2.,  0.,  2.]]]], shape=(4392, 2, 16, 8))

datasets = remove_incomplete_days(dat,T)
datasets

(array([[[[ 0.,  0.,  1., ...,  0.,  0.,  0.],
          [ 0.,  1.,  1., ...,  1.,  0.,  0.],
          [ 0.,  0.,  2., ...,  2.,  2.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  3.,  1.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  1.]],
 
         [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  0.,  1.,  0.],
          [ 0.,  0.,  0., ...,  2.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.]]],
 
 
        [[[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  0.,  1.,  2.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  1.]],
 
         [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  1., ...,  0.,  0.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  1.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  1.]]],
 
 
        [[[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  0.,  2.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.]],
 
         [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          [ 0.,  0.,  0., ...,  2.,  1.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.]]],
 
 
        ...,
 
 
        [[[ 0.,  0.,  5., ...,  0.,  0.,  0.],
          [ 0.,  2.,  6., ...,  6.,  0.,  0.],
          [ 0.,  5., 22., ..., 14.,  8.,  0.],
          ...,
          [ 0.,  0.,  4., ...,  9.,  5.,  2.],
          [ 0.,  0.,  3., ...,  7.,  3.,  1.],
          [ 0.,  0.,  0., ...,  4.,  6.,  8.]],
 
         [[ 0.,  0.,  2., ...,  0.,  0.,  0.],
          [ 0.,  2.,  2., ...,  9.,  0.,  0.],
          [ 0.,  0., 10., ..., 13.,  7.,  0.],
          ...,
          [ 0.,  0.,  3., ...,  7.,  4.,  1.],
          [ 0.,  0.,  2., ...,  3.,  2.,  2.],
          [ 0.,  0.,  0., ...,  2.,  1.,  1.]]],
 
 
        [[[ 0.,  0.,  5., ...,  0.,  0.,  0.],
          [ 0.,  1.,  4., ...,  5.,  0.,  0.],
          [ 0.,  1.,  7., ...,  3.,  7.,  0.],
          ...,
          [ 0.,  0.,  3., ...,  8.,  1.,  0.],
          [ 0.,  0.,  2., ...,  4.,  0.,  1.],
          [ 0.,  0.,  0., ...,  2.,  4.,  3.]],
 
         [[ 0.,  0.,  2., ...,  0.,  0.,  0.],
          [ 0.,  1.,  5., ...,  7.,  0.,  0.],
          [ 0.,  0.,  8., ...,  3.,  5.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  2.,  1.,  1.],
          [ 0.,  0.,  2., ...,  0.,  1.,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  1.]]],
 
 
        [[[ 0.,  0.,  1., ...,  0.,  0.,  0.],
          [ 0.,  0.,  3., ...,  4.,  0.,  0.],
          [ 0.,  2.,  2., ...,  4.,  6.,  0.],
          ...,
          [ 0.,  0.,  0., ..., 12.,  0.,  2.],
          [ 0.,  0.,  1., ...,  3.,  1.,  0.],
          [ 0.,  0.,  0., ...,  2.,  2.,  2.]],
 
         [[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  2., ...,  3.,  0.,  0.],
          [ 0.,  0.,  5., ...,  0.,  2.,  0.],
          ...,
          [ 0.,  0.,  0., ...,  1.,  0.,  0.],
          [ 0.,  0.,  1., ...,  1.,  0.,  0.],
          [ 0.,  0.,  0., ...,  2.,  0.,  2.]]]], shape=(4392, 2, 16, 8)),
 array(['2014040101', '2014040102', '2014040103', ..., '2014093022',
        '2014093023', '2014093024'], shape=(4392,), dtype='<U10'))

datasets[1]

array(['2014040101', '2014040102', '2014040103', ..., '2014093022',
       '2014093023', '2014093024'], shape=(4392,), dtype='<U10')

def extract_timeseries(data: np.ndarray, row: int, col: int, channel: int = 0) -> np.ndarray:
    """
    Extracts the time series of a specific grid location from a 4D numpy array.
    
    Parameters:
        data (np.ndarray): A 4D numpy array with shape (time_steps, channels, height, width)
        row (int): The row index of the desired grid location
        col (int): The column index of the desired grid location
        channel (int): The channel index to extract from (default: 0)
        
    Returns:
        np.ndarray: A 1D array representing the time series at the specified location
    """
    if data.ndim != 4:
        raise ValueError("Input data must be a 4D numpy array (time_steps, channels, height, width)")
    
    if not (0 <= row < data.shape[2]) or not (0 <= col < data.shape[3]):
        raise ValueError("Row or column index out of bounds")
    
    if not (0 <= channel < data.shape[1]):
        raise ValueError("Channel index out of bounds")
    
    return data[:, channel, row, col]

grid_data = datasets[0]

ROW=3
COL=5

# Extract time series from row 5, column 3 in channel 0
time_series = extract_timeseries(grid_data, row=ROW, col=COL, channel=0)

print(time_series)

[0. 1. 0. ... 2. 4. 1.]

len(time_series)

4392

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

def adjust_timestamps(timestamps):
    """
    Adjusts the given timestamps by subtracting 1 hour from each entry.
    
    Parameters:
        timestamps (numpy array): Array of timestamps in 'YYYYMMDDHH' format.
        
    Returns:
        numpy array: Adjusted timestamps.
    """
    adjusted = []
    for ts in timestamps:
        # Fix invalid '24' hour case
        if ts[-2:] == '24':
            ts = (datetime.strptime(ts[:8], "%Y%m%d") + timedelta(days=1)).strftime("%Y%m%d") + "00"
        dt = datetime.strptime(ts, "%Y%m%d%H") - timedelta(hours=1)
        adjusted.append(dt.strftime("%Y%m%d%H"))
    
    return np.array(adjusted, dtype='<U10')

# Example usage
timestamps = datasets[1]
timestamps = adjust_timestamps(timestamps)
timestamps

array(['2014040100', '2014040101', '2014040102', ..., '2014093021',
       '2014093022', '2014093023'], shape=(4392,), dtype='<U10')

def format_timestamp(timestamp: str) -> str:
    """
    Converts timestamp from 'YYYYMMDDHH' to 'YYYY-MM-DD HH:00'.
    """
    dt = datetime.strptime(timestamp, "%Y%m%d%H")
    return dt.strftime("%Y-%m-%d %H:00")

def create_dataframe(timeseries: np.ndarray, timestamps: np.ndarray) -> pd.DataFrame:
    if timeseries.shape[0] != timestamps.shape[0]:
        raise ValueError("Time series and timestamps must have the same length")

    formatted_timestamps = [format_timestamp(ts) for ts in timestamps]

    df = pd.DataFrame({
        'timestamp': formatted_timestamps,
        'value': timeseries
    })
    return df


df = create_dataframe(time_series, timestamps)
print(df)

             timestamp  value
0     2014-04-01 00:00    0.0
1     2014-04-01 01:00    1.0
2     2014-04-01 02:00    0.0
3     2014-04-01 03:00    0.0
4     2014-04-01 04:00    1.0
...                ...    ...
4387  2014-09-30 19:00   16.0
4388  2014-09-30 20:00   10.0
4389  2014-09-30 21:00    2.0
4390  2014-09-30 22:00    4.0
4391  2014-09-30 23:00    1.0

[4392 rows x 2 columns]

train_data_df = df[df['timestamp']<'2014-09-21 00:00']

test_data_df = df.copy()

df = df[df['timestamp']>='2014-09-16 00:00']

from nixtla import NixtlaClient

nixtla_client = NixtlaClient(
    api_key = "nixak-"
)

nixtla_client.validate_api_key()

INFO:nixtla.nixtla_client:Happy Forecasting! :)

True

nixtla_client.plot(df, time_col='timestamp', target_col='value')

timegpt_fcst_df = nixtla_client.forecast(df=df, h=24, freq='h', time_col='timestamp', target_col='value')
timegpt_fcst_df

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

nixtla_client.plot(df, timegpt_fcst_df, time_col='timestamp', target_col='value')

timegpt_fcst_with_history_df = nixtla_client.forecast(
    df=df, h=1, time_col='timestamp', target_col='value', freq='h',
    add_history=True,
)

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...
INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...

timegpt_fcst_with_history_df.to_csv(f'predictions/timegpt/fcst_{ROW}_{COL}_h{HORIZON}.csv')

nixtla_client.plot(df, timegpt_fcst_with_history_df, time_col='timestamp', target_col='value')

def compute_rmse(original_df, forecast_df):
    # Convert timestamp to datetime
    original_df['timestamp'] = pd.to_datetime(original_df['timestamp'])
    forecast_df['timestamp'] = pd.to_datetime(forecast_df['timestamp'])

    # Merge the DataFrames based on the timestamp
    merged_df = pd.merge(original_df, forecast_df, on='timestamp', how='inner')
    print(merged_df)

    # Calculate the squared differences
    merged_df['squared_diff'] = (merged_df['value'] - merged_df['TimeGPT'])**2

    # Calculate RMSE
    rmse = np.sqrt(merged_df['squared_diff'].mean())
    return rmse

rmse = compute_rmse(df,timegpt_fcst_with_history_df)
print(f"RMSE: {rmse}")

import os.path

def prediction_exists(row, col):
    return os.path.isfile(f'predictions/timegpt/fcst_{row}_{col}.csv')

prediction_exists(2,1)

False

def load_prediction(row,col):
    timegpt_fcst_with_history_df = pd.read_csv(f'predictions/timegpt/fcst_{row}_{col}.csv')
    return timegpt_fcst_with_history_df

load_prediction(1,3)

def compute_cell_rmse(row, col):
    time_series = extract_timeseries(grid_data, row=row, col=col, channel=0)
    timestamps = datasets[1]
    timestamps = adjust_timestamps(timestamps)
    df = create_dataframe(time_series, timestamps)
    df = df[df['timestamp']>='2014-09-16 00:00']

    if prediction_exists(row, col):
        print('Loading predictions ...')
        timegpt_fcst_with_history_df = load_prediction(row, col)
    else:
        print('Computing predicitons ...')
        timegpt_fcst_with_history_df = nixtla_client.forecast(
            df=df, h=1, time_col='timestamp', target_col='value', freq='h',
            add_history=True,
        )
        timegpt_fcst_with_history_df.to_csv(f'predictions/timegpt/fcst_{row}_{col}.csv', index=False)
    rmse = compute_rmse(df, timegpt_fcst_with_history_df)
    print(f"[{row},{col}] RMSE:{rmse}")
    return rmse

rmses =  []

for row in range(0,16):
    for col in range(0,8):
        rmses.append(compute_cell_rmse(row,col))

print(rmses)

Loading predictions ...
[0,0] RMSE:0.002570502930800306
Loading predictions ...
[0,1] RMSE:0.002570502930800306
Loading predictions ...
[0,2] RMSE:3.5633947837540783
Loading predictions ...
[0,3] RMSE:10.037934923310697
Loading predictions ...
[0,4] RMSE:15.239814706642761
Loading predictions ...
[0,5] RMSE:0.002570502930800306

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

Loading predictions ...
[0,6] RMSE:0.002570502930800306
Loading predictions ...
[0,7] RMSE:0.002570502930800306
Computing predicitons ...

INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...
INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

[1,0] RMSE:0.002570502930800306
Loading predictions ...
[1,1] RMSE:6.072188987583208
Computing predicitons ...

INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...
INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

[1,2] RMSE:6.10467960056596
Loading predictions ...
[1,3] RMSE:10.603164877491775
Computing predicitons ...

INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...
INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

[1,4] RMSE:13.917665751860108
Computing predicitons ...

INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...
INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

np.mean(rmses)

np.float64(5.7716825988548335)

def get_df(row, col, i=0, channel=0):
    time_series = extract_timeseries(grid_data, row=row, col=col, channel=channel)
    timestamps = datasets[1]
    timestamps = adjust_timestamps(timestamps)
    df = create_dataframe(time_series, timestamps)
    df['unique_id'] = i
    return df

dfs = []
i=0

for channel in [0,1]:
    for row, col in [(3,5),(4,5)]:
        dfs.append(get_df(row,col,i,channel))
        i+=1

df = pd.concat(dfs, axis=0)
df

nixtla_client.plot(df, time_col='timestamp', target_col='value')

timegpt_fcst_df = nixtla_client.forecast(df=df, h=24, freq='h', time_col='timestamp', target_col='value')
timegpt_fcst_df

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Querying model metadata...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

nixtla_client.plot(df, timegpt_fcst_df, time_col='timestamp', target_col='value', max_insample_length=24*11)

def backtest_forecasting(df, start_timestamp, horizon, timestamp='ds', target_col='y'):
    """
    Backtest a time forecasting model by iteratively making forecasts at fixed intervals.
    
    :param df: DataFrame with time series data (must contain timestamp and target_col columns)
    :param start_timestamp: The timestamp from which to start the backtesting process
    :param horizon: Forecasting horizon (number of time steps to forecast at each iteration)
    
    :return: A DataFrame with concatenated forecasts
    """
    
    # Ensure timestamp is in datetime format
    df[timestamp] = pd.to_datetime(df[timestamp])
    
    all_forecasts = []
    current_timestamp = start_timestamp

    print(f'current: {current_timestamp}')
    print(f'max: {df[timestamp].max()}')
    
    while current_timestamp < df[timestamp].max():
        # Filter the data up to the current timestamp

        print(f'current: {current_timestamp}')
        df_subset = df[df[timestamp] < current_timestamp]
        
        # Make a forecast using the model
        forecast_df = nixtla_client.forecast(df=df_subset, h=horizon, freq='h', time_col=timestamp, target_col=target_col)
        
        # Append the results
        all_forecasts.append(forecast_df)
        
        # Move to the next forecasting window
        current_timestamp += pd.Timedelta(hours=horizon)
    
    # Concatenate all forecasts into a single DataFrame
    return pd.concat(all_forecasts, ignore_index=True)

timegpt_fcst_df_backtest = backtest_forecasting(df, datetime(2014,9,21), 24, timestamp='timestamp', target_col='value')

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-21 00:00:00
max: 2014-09-30 23:00:00
current: 2014-09-21 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-22 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-23 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-24 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-25 00:00:00

nixtla_client.plot(df, timegpt_fcst_df_backtest, time_col='timestamp', target_col='value', max_insample_length=24*11)

rmse = compute_rmse(df, timegpt_fcst_df_backtest)
rmse

               timestamp  value  unique_id_x  unique_id_y   TimeGPT
0    2014-09-21 01:00:00    4.0            0            0  2.516364
1    2014-09-21 01:00:00    4.0            0            1 -0.000597
2    2014-09-21 01:00:00    4.0            0            2  0.754800
3    2014-09-21 01:00:00    4.0            0            3 -1.170985
4    2014-09-21 02:00:00    1.0            0            0  1.462680
...                  ...    ...          ...          ...       ...
3819 2014-09-30 22:00:00    7.0            3            3  2.302446
3820 2014-09-30 23:00:00    5.0            3            0  5.505838
3821 2014-09-30 23:00:00    5.0            3            1  4.514214
3822 2014-09-30 23:00:00    5.0            3            2  3.481431
3823 2014-09-30 23:00:00    5.0            3            3  2.236673

[3824 rows x 5 columns]

np.float64(8.346675647024847)

timegpt_fcst_df_backtest = backtest_forecasting(df, datetime(2014,9,21), 48, timestamp='timestamp', target_col='value')

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
WARNING:nixtla.nixtla_client:The specified horizon "h" exceeds the model horizon, this may lead to less accurate forecasts. Please consider using a smaller horizon.
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-21 00:00:00
max: 2014-09-30 23:00:00
current: 2014-09-21 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
WARNING:nixtla.nixtla_client:The specified horizon "h" exceeds the model horizon, this may lead to less accurate forecasts. Please consider using a smaller horizon.
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-23 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
WARNING:nixtla.nixtla_client:The specified horizon "h" exceeds the model horizon, this may lead to less accurate forecasts. Please consider using a smaller horizon.
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-25 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
WARNING:nixtla.nixtla_client:The specified horizon "h" exceeds the model horizon, this may lead to less accurate forecasts. Please consider using a smaller horizon.
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-27 00:00:00

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
WARNING:nixtla.nixtla_client:The specified horizon "h" exceeds the model horizon, this may lead to less accurate forecasts. Please consider using a smaller horizon.
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

current: 2014-09-29 00:00:00

nixtla_client.plot(df, timegpt_fcst_df_backtest, time_col='timestamp', target_col='value', max_insample_length=24*11)

rmse = compute_rmse(df, timegpt_fcst_df_backtest)
rmse

               timestamp  value  unique_id_x  unique_id_y   TimeGPT
0    2014-09-21 01:00:00    4.0            0            0  2.516364
1    2014-09-21 01:00:00    4.0            0            1 -0.000597
2    2014-09-21 01:00:00    4.0            0            2  0.754800
3    2014-09-21 01:00:00    4.0            0            3 -1.170985
4    2014-09-21 02:00:00    1.0            0            0  1.462680
...                  ...    ...          ...          ...       ...
3819 2014-09-30 22:00:00    7.0            3            3  2.953221
3820 2014-09-30 23:00:00    5.0            3            0  3.396904
3821 2014-09-30 23:00:00    5.0            3            1  4.583953
3822 2014-09-30 23:00:00    5.0            3            2  2.344816
3823 2014-09-30 23:00:00    5.0            3            3  2.392262

[3824 rows x 5 columns]

np.float64(8.666202218206987)

def compute_rmse(original_df, forecast_df, timestamp='ds', target_col='y', predict_col='TimeGPT', unique_id='unique_id'):
    # Convert timestamp to datetime
    original_df[timestamp] = pd.to_datetime(original_df[timestamp])
    forecast_df[timestamp] = pd.to_datetime(forecast_df[timestamp])

    # Merge the DataFrames based on the timestamp
    merged_df = pd.merge(original_df, forecast_df, on=[timestamp,unique_id], how='inner')
    #print(merged_df)

    # Calculate the squared differences
    merged_df['squared_diff'] = (merged_df[target_col] - merged_df[predict_col])**2

    # Calculate RMSE
    rmse = np.sqrt(merged_df['squared_diff'].mean())
    return rmse

dfs = []
i=0

for channel in [0,1]:
    for row in range(0,16):
        for col in range(0,8):
            dfs.append(get_df(row,col,i,channel))
            i+=1

df = pd.concat(dfs, axis=0)
df

import os 

HORIZON = 1
path = f'predictions/timegpt/fcst_c01_h{HORIZON}.csv'

if os.path.isfile(path):
    timegpt_fcst_df_backtest = pd.read_csv(path)
else: 
    timegpt_fcst_df_backtest = backtest_forecasting(df, datetime(2014,9,21), HORIZON, timestamp='timestamp', target_col='value')
    timegpt_fcst_df_backtest.to_csv(path, index=False)

rmse = compute_rmse(df, timegpt_fcst_df_backtest, timestamp='timestamp', target_col='value')
print(f'horizon={HORIZON}: rmse={rmse}')

                timestamp  value  unique_id   TimeGPT
0     2014-09-21 01:00:00    0.0          0  0.003328
1     2014-09-21 02:00:00    0.0          0  0.003328
2     2014-09-21 03:00:00    0.0          0  0.003328
3     2014-09-21 04:00:00    0.0          0  0.003328
4     2014-09-21 05:00:00    0.0          0  0.003328
...                   ...    ...        ...       ...
61179 2014-09-30 19:00:00    1.0        255  3.210955
61180 2014-09-30 20:00:00    4.0        255  1.871202
61181 2014-09-30 21:00:00    1.0        255  2.555601
61182 2014-09-30 22:00:00    1.0        255  0.832277
61183 2014-09-30 23:00:00    2.0        255  0.272251

[61184 rows x 4 columns]
horizon=1: rmse=5.695230340297138

nixtla_client.plot(df, timegpt_fcst_df_backtest, time_col='timestamp', target_col='value', max_insample_length=24*11, unique_ids=[22,36,51,60,86,96,108,127])

import pandas as pd
from nixtla import NixtlaClient

nixtla_client = NixtlaClient(
    api_key = "nixak-"
)
nixtla_client.validate_api_key()

INFO:nixtla.nixtla_client:Happy Forecasting! :)

True

df = pd.read_csv('datasets/BikeVienna/service_2-2019-05-07_2019-12-31.csv', usecols=['timestamp','station_id','vehicles_available'])
df.rename(columns={'timestamp':'ds','station_id':'unique_id', 'vehicles_available':'y'}, inplace=True)
df['ds'] = pd.to_datetime(df['ds']).dt.tz_convert(None)
df

df = df[df.ds<'2019-08-15']
df

len(df.unique_id.unique())

120

nixtla_client.plot(df, max_insample_length=2400)

df = df.set_index('ds').groupby('unique_id').resample('60min')['y'].max()
df.ffill(inplace=True)
df = df.reset_index()
df

nixtla_client.plot(df, max_insample_length=2400)

from statsforecast import StatsForecast
from statsforecast.models import (
    HistoricAverage,
    SeasonalNaive,
    AutoARIMA
)

models = [
    HistoricAverage(),
    SeasonalNaive(season_length=24),
    AutoARIMA(),
    
]

sf = StatsForecast( 
    models=models,
    freq='h', 
    fallback_model = SeasonalNaive(season_length=7),
    n_jobs=-1,
)

forecasts_df = sf.forecast(df=df, h=24, level=[])
forecasts_df

nixtla_client.plot(df, forecasts_df, max_insample_length=24*11)

def backtest_baselines(df, start_timestamp, horizon, timestamp='ds', target_col='y'):
    df[timestamp] = pd.to_datetime(df[timestamp])
    
    all_forecasts = []
    current_timestamp = start_timestamp

    print(f'current: {current_timestamp}')
    print(f'max: {df[timestamp].max()}')
    
    while current_timestamp < df[timestamp].max():
        print(f'current: {current_timestamp}')
        df_subset = df[df[timestamp] <= current_timestamp]
        
        forecast_df = sf.forecast(df=df_subset, h=horizon, level=[], time_col=timestamp, target_col=target_col)
        all_forecasts.append(forecast_df)
        
        current_timestamp += pd.Timedelta(hours=horizon)
    
    return pd.concat(all_forecasts, ignore_index=True)

import os
import numpy as np
from datetime import datetime

HORIZON = 12
path = f'predictions/timegpt/baselines_vienna_h{HORIZON}.csv'

if os.path.isfile(path):
    baseline_df_backtest = pd.read_csv(path)
else: 
    baseline_df_backtest = backtest_baselines(df, datetime(2019,8,5), HORIZON)
    baseline_df_backtest.to_csv(path, index=False)

for model_name in ['HistoricAverage','SeasonalNaive','AutoARIMA']:
    rmse = compute_rmse(df, baseline_df_backtest, predict_col=model_name)		
    print(f'{model_name} horizon={HORIZON}: rmse={rmse}')

HistoricAverage horizon=12: rmse=7.035121951753955
SeasonalNaive horizon=12: rmse=7.870333843727055
AutoARIMA horizon=12: rmse=5.2139226565300625

nixtla_client.plot(df, baseline_df_backtest, max_insample_length=24*11)

tmp = df[df.unique_id==16]

timegpt_fcst_df = nixtla_client.forecast(df=tmp, h=24, freq='h')
timegpt_fcst_df

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Restricting input...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...

nixtla_client.plot(tmp, timegpt_fcst_df,  max_insample_length=24*15)

timegpt_fcst_with_history_df = nixtla_client.forecast(
    df=tmp, h=1, freq='h', add_history=True,
)

INFO:nixtla.nixtla_client:Validating inputs...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...
INFO:nixtla.nixtla_client:Preprocessing dataframes...
INFO:nixtla.nixtla_client:Calling Forecast Endpoint...
INFO:nixtla.nixtla_client:Calling Historical Forecast Endpoint...

nixtla_client.plot(tmp, timegpt_fcst_with_history_df)

df = df[df.ds<'2019-08-15']
df

import os 
import numpy as np
from datetime import datetime

HORIZON = 12
path = f'predictions/timegpt/fcst_vienna_h{HORIZON}.csv'

if os.path.isfile(path):
    timegpt_fcst_df_backtest = pd.read_csv(path)
else: 
    timegpt_fcst_df_backtest = backtest_forecasting(df, datetime(2019,8,5), HORIZON)
    timegpt_fcst_df_backtest.to_csv(path, index=False)

rmse = compute_rmse(df, timegpt_fcst_df_backtest)
print(f'horizon={HORIZON}: rmse={rmse}')

horizon=12: rmse=4.838292942837482

nixtla_client.plot(df, timegpt_fcst_df_backtest, max_insample_length=24*11)

	timestamp	TimeGPT
0	2014-10-01 00:00:00	-1.202475
1	2014-10-01 01:00:00	-0.848185
2	2014-10-01 02:00:00	-0.593907
3	2014-10-01 03:00:00	-1.358113
4	2014-10-01 04:00:00	-2.240385
5	2014-10-01 05:00:00	-1.024400
6	2014-10-01 06:00:00	9.498710
7	2014-10-01 07:00:00	28.237295
8	2014-10-01 08:00:00	53.401657
9	2014-10-01 09:00:00	40.392822
10	2014-10-01 10:00:00	16.469150
11	2014-10-01 11:00:00	9.792094
12	2014-10-01 12:00:00	12.540842
13	2014-10-01 13:00:00	13.419383
14	2014-10-01 14:00:00	12.166123
15	2014-10-01 15:00:00	10.900003
16	2014-10-01 16:00:00	12.008976
17	2014-10-01 17:00:00	17.562895
18	2014-10-01 18:00:00	18.652906
19	2014-10-01 19:00:00	14.093687
20	2014-10-01 20:00:00	8.132493
21	2014-10-01 21:00:00	4.112089
22	2014-10-01 22:00:00	3.221751
23	2014-10-01 23:00:00	2.066868

	unique_id	timestamp	TimeGPT
0	0	2014-10-01 00:00:00	-1.202614
1	0	2014-10-01 01:00:00	-0.847775
2	0	2014-10-01 02:00:00	-0.594713
3	0	2014-10-01 03:00:00	-1.358931
4	0	2014-10-01 04:00:00	-2.239870
...	...	...	...
91	3	2014-10-01 19:00:00	25.792120
92	3	2014-10-01 20:00:00	12.507367
93	3	2014-10-01 21:00:00	8.102301
94	3	2014-10-01 22:00:00	6.962522
95	3	2014-10-01 23:00:00	5.438268

Model	RMSE
ARIMA	10.56
SARIMA	10.07
VAR	9.92
DeepST-C	8.39
DeepST-CP	7.64
DeepST-CPT	7.56
DeepST-CPTM	7.43
ST-ResNet	6.33
TimeGPT (horizon=1)	5.70
TimeGPT (horizon=12)	7.62
TimeGPT (horizon=24)	8.93

	ds	unique_id	y
0	2019-05-06 22:07:01	53	4
1	2019-05-06 22:24:01	53	4
2	2019-05-06 22:37:01	53	4
3	2019-05-06 22:52:01	53	4
4	2019-05-06 23:07:01	53	4
...	...	...	...
2538844	2019-12-31 21:52:01	128	22
2538845	2019-12-31 22:07:01	128	22
2538846	2019-12-31 22:24:01	128	22
2538847	2019-12-31 22:37:01	128	22
2538848	2019-12-31 22:52:01	128	22

	ds	unique_id	y
0	2019-05-06 22:07:01	53	4
1	2019-05-06 22:24:01	53	4
2	2019-05-06 22:37:01	53	4
3	2019-05-06 22:52:01	53	4
4	2019-05-06 23:07:01	53	4
...	...	...	...
2527148	2019-08-14 22:52:01	128	6
2527149	2019-08-14 23:07:01	128	11
2527150	2019-08-14 23:24:01	128	11
2527151	2019-08-14 23:37:01	128	11
2527152	2019-08-14 23:52:01	128	11

The Future of Urban Mobility AI: Can Foundation Models Outperform Classic ML & DL Models For Crowd Prediction?¶

Original Paper: Zhang et al. (2017) "Deep Spatio-Temporal Residual Networks for Citywide Crowd Flows Prediction" (Microsoft Research)¶

Goal¶

Datasets¶

Results¶

Code¶

BikeNYC Experiment¶

TimeGPT¶

Experiment v0¶

Simple forecast¶

Historical forecast¶

Run forecasts for all grid cells and store results¶

Identified issues¶

Experiment v1¶

Multiple cells at once¶

Backtesting¶

All cells & backtesting¶

Results¶

Conclusion¶

BikeVienna Experiment¶

Resampling¶

Baseline¶

TimeGPT¶

Simple forecast¶

Historical forecast¶

All stations and backtesting¶

Results¶

Remarks¶

Earlier paper by Zhang et al. (2016) "DNN-Based Prediction Model for Spatio-Temporal Data"¶

Multi-step forcasting (i.e. longer forecasting horizon)¶

More data is better ;-)¶

Potential error in RMSE reporting¶

	timestamp	TimeGPT
0	2014-09-21 00:00:00	16.648037
1	2014-09-21 01:00:00	10.752392
2	2014-09-21 02:00:00	6.054309
3	2014-09-21 03:00:00	2.488530
4	2014-09-21 04:00:00	1.055259
...	...	...
236	2014-09-30 20:00:00	28.864220
237	2014-09-30 21:00:00	26.119790
238	2014-09-30 22:00:00	21.338356
239	2014-09-30 23:00:00	13.381630
240	2014-10-01 00:00:00	9.842307

	timestamp	value	unique_id
0	2014-04-01 00:00	0.0	0
1	2014-04-01 01:00	1.0	0
2	2014-04-01 02:00	0.0	0
3	2014-04-01 03:00	0.0	0
4	2014-04-01 04:00	1.0	0
...	...	...	...
4387	2014-09-30 19:00	29.0	3
4388	2014-09-30 20:00	8.0	3
4389	2014-09-30 21:00	6.0	3
4390	2014-09-30 22:00	7.0	3
4391	2014-09-30 23:00	5.0	3

	unique_id	ds	y
0	9	2019-05-06 22:00:00	8.0
1	9	2019-05-06 23:00:00	8.0
2	9	2019-05-07 00:00:00	8.0
3	9	2019-05-07 01:00:00	8.0
4	9	2019-05-07 02:00:00	8.0
...	...	...	...
288235	128	2019-08-14 19:00:00	5.0
288236	128	2019-08-14 20:00:00	4.0
288237	128	2019-08-14 21:00:00	4.0
288238	128	2019-08-14 22:00:00	6.0
288239	128	2019-08-14 23:00:00	11.0

	unique_id	ds	HistoricAverage	SeasonalNaive	AutoARIMA
0	9	2019-08-15 00:00:00	14.399667	7.0	17.722890
1	9	2019-08-15 01:00:00	14.399667	5.0	17.505003
2	9	2019-08-15 02:00:00	14.399667	5.0	17.245954
3	9	2019-08-15 03:00:00	14.399667	6.0	17.099813
4	9	2019-08-15 04:00:00	14.399667	6.0	17.014962
...	...	...	...	...	...
2875	128	2019-08-15 19:00:00	15.333472	5.0	9.244525
2876	128	2019-08-15 20:00:00	15.333472	4.0	9.435932
2877	128	2019-08-15 21:00:00	15.333472	4.0	9.637041
2878	128	2019-08-15 22:00:00	15.333472	6.0	9.782622
2879	128	2019-08-15 23:00:00	15.333472	11.0	9.941062

	unique_id	ds	TimeGPT
0	16	2019-12-31 23:00:00	17.982372
1	16	2020-01-01 00:00:00	17.792389
2	16	2020-01-01 01:00:00	17.701164
3	16	2020-01-01 02:00:00	17.599297
4	16	2020-01-01 03:00:00	17.448940
5	16	2020-01-01 04:00:00	17.369642
6	16	2020-01-01 05:00:00	17.177364
7	16	2020-01-01 06:00:00	16.986202
8	16	2020-01-01 07:00:00	16.813074
9	16	2020-01-01 08:00:00	16.708633
10	16	2020-01-01 09:00:00	16.464092
11	16	2020-01-01 10:00:00	16.334290
12	16	2020-01-01 11:00:00	16.334015
13	16	2020-01-01 12:00:00	16.260014
14	16	2020-01-01 13:00:00	16.124367
15	16	2020-01-01 14:00:00	15.809680
16	16	2020-01-01 15:00:00	15.547085
17	16	2020-01-01 16:00:00	15.324635
18	16	2020-01-01 17:00:00	15.201200
19	16	2020-01-01 18:00:00	15.067950
20	16	2020-01-01 19:00:00	14.818552
21	16	2020-01-01 20:00:00	14.678573
22	16	2020-01-01 21:00:00	14.522839
23	16	2020-01-01 22:00:00	14.376727

Model	Horizon	RMSE
HistoricAverage	1	7.0229
HistoricAverage	12	7.0195
HistoricAverage	24	7.0426
SeasonalNaive	1	7.8703
SeasonalNaive	12	7.7317
SeasonalNaive	24	7.8703
AutoARIMA	1	2.2639
AutoARIMA	12	5.1505
AutoARIMA	24	6.3881
TimeGPT	1	2.3193
TimeGPT	12	4.8383
TimeGPT	24	5.6671